Transcriptomic profiling and discovery of key transcription factors involved in adventitious roots formation from root cuttings of mulberry

ARs plays a crucial role in plant morphogenesis and development. The limited and inefficient rooting of scions poses a significant challenge to the efficiency and quality of clonal propagation of forest trees in silvicultural practices. Building on previous research conducted by our team, we found that applying IBA at a concentration of 1000 mg/L significantly enhanced mulberry rooting. This study aims to uncover the molecular mechanisms underlying this effect by analyzing RNA sequencing data from mulberry phloem before and after treatment with IBA over time intervals of 10, 20, 30, and 40 days. We identified 5226 DEGs, which were then classified into GO terms and KEGG pathways, showing significant enrichment in hormone signaling processes. Using WGCNA, we identified eight co-expression modules, two of which were significantly correlated with the IBA treatment. Additionally, 18 transcription factors that potentially facilitate ARs formation in mulberry were identified, and an exploratory analysis on the cis-regulatory elements associated with these transcription factors was conducted. The findings of this study provide a comprehensive understanding of the mechanisms of ARs in mulberry and offer theoretical support for the discovery and utilization of exceptional genetic resources within the species.


Introduction
Adventitious roots (ARs), formed from non-root organs such as stems and leaves [1], enhance a plant's ability to adapt to environmental changes and play a vital role in plant morphogenesis and development [2,3].The development of ARs in woody plants can be divided into three stages: dedifferentiation, induction, and differentiation scarcity of root formation are pivotal issues that impact the efficiency and quality of clonal propagation of forest trees in forestry operations [8].
Numerous studies have demonstrated that various phytohormones significantly influence the formation and development of ARs, with auxins, particularly indole-3-acetic acid (IAA) [9], having the most profound effect.Pei Dong and colleagues [10] have proposed that elevated IAA levels play a crucial role in the differentiation of root primordia during the initial phases of ARs induction.During this phase, IAA concentrations increase at the scion incision site, but these levels decrease once the root primordia are established [11].Besides IAA, other phytohormones such as ethylene also play a role in AR formation, as demonstrated by its ability to stimulate root regeneration in species including chrysanthemum, petunia, and Arabidopsis thaliana [12][13][14].The interaction between ethylene and IAA can synergistically enhance ARs genesis [15][16][17].It has been shown that genes from the ethylene-responsive AP2/ERF transcription factor family are upregulated during the ARs induction period in poplar [18].Conversely, cytokinins may antagonize auxin activity and inhibit AR development across various plant species [19][20][21], with higher IAA/cytokinin ratios being conducive to AR formation.Wang and collaborators [22] have found that abscisic acid can promote rooting in tetraploid acacia scions by counteracting the suppressive effects of high IAA concentrations.Additionally, Gutierrez and team [23] have shown that the auxin-responsive Gretchen Hagen3 (GH3) gene family, specifically GH3.3, GH3.5, and GH3.6, are crucial for the fine-tuning of ARs initiation in Arabidopsis through the modulation of jasmonic acid homeostasis.
Advances in ARs research have moved from anatomical and physiological studies to the molecular level, largely driven by the development and integration of RNA sequencing (RNA-seq) technology.RNA-seq analysis has revealed that phytohormone signaling pathways are predominant in ARs development, as indicated in studies by Li Ke [24], who noted that exogenous indole-3-butyric acid (IBA) significantly induced hormone biosynthesis and responsive gene expression during ARs development in apple rootstocks.Similarly, Cheng Long's RNA-seq studies [25] suggested that aluminum exposure might facilitate the regeneration and development of ARs in tea plants through a complex transcriptional regulatory network involving various plant hormones and associated genes.However, research on the molecular mechanisms of ARs formation has primarily focused on model plants like Arabidopsis and rice, with limited studies on mulberry.
Mulberry has been cultivated in China for a long time, and with advancements in research methodologies, our understanding of its properties has deepened.Mulberry leaves and fruits are known for their high nutritional value and health benefits [26][27][28].Although mulberry cuttings traditionally exhibit low survival rates, advancements in cutting techniques have significantly improved their viability [29,30].However, the generation and quality of roots remain major challenges.Quickly forming mature roots in cutting seedlings is currently an urgent issue to address [31].
According to prior research conducted by our group, the application of IBA at a concentration of 1000 mg/L was found to be most effective for mulberry rooting [32].To explore the molecular mechanisms, RNA-seq was performed at intervals before and after the 1000 mg/L IBA treatment, leading to the identification of differentially expressed genes (DEGs) categorized into Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.Subsequently, a transcription factor (TF) gene regulatory network was constructed from these DEGs using Weighted Gene Correlation Network Analysis (WGCNA).The enrichment of DEGs in GO terms and KEGG pathways, along with the construction of the TF gene regulatory network based on WGCNA, aims to provide both practical and theoretical insights for the propagation and rooting mechanisms of mulberry cuttings, potentially benefiting the cultivation practices of other woody plant species.

Experimental site, plant materials and experimental design
The experiment was carried out at the third living area of Henan Agricultural University located in Zhengzhou City, Henan Province.The site is situated at geographical coordinates of 113.22°E longitude, 34.28°N latitude, with an elevation of 98 m above sea level.The cuttings were sourced from a greenhouse specifically designed for propagation, which was equipped with comprehensive full-spectrum lighting and an automatic misting system.The greenhouse contained cutting pools divided into five sections, each approximately 11 m by 6 m, with each subpool measuring 6 m by 2 m and having a depth of 0.4 m.
Cuttings were taken from semi-lignified branches of the mulberry cultivar "Qiangsang No. 1, " developed by the Silkworm Research Institute of the Zhejiang Academy of Agricultural Sciences.These were processed into uniformly sized stakes, and the basal ends were treated with either a 1000 mg/L IBA solution for the treatment group or water for the control group (CK) for 30 s.Each treatment was replicated three times, using 80 cuttings per replicate, all cuttings were guaranteed to come from a uniform clone Soaked all the branches in carbendazim for 1 ~ 2 min, put the cuttings in a cool and ventilated place, dry the liquid in the shade.The prepared cuttings were then inserted into a sterilized growth medium according to established protocols detailed in our previous publication [32].The developmental stages of the softwood cuttings' rooting process were documented, as illustrated in Fig. 1.This documentation included stages of callus formation from 0 to 10 days post-planting, induction of root primordia from 10 to 20 days, emergence and formation of ARs from 20 to 30 days, and elongation and maturation of ARs beyond 40 days.For the study of ARs development, cortical tissue samples approximately 1 cm above the base of the cuttings were harvested at 10 (CK-1, IBA-1), 20 (CK-2, IBA-2), 30 (CK-3, IBA-3), and 40 days (CK-4, IBA-4) post-planting for both the control and treatment groups.The collected samples were immediately immersed in liquid nitrogen and subsequently stored in a -80 °C freezer.Twenty specimens were randomly selected from each time point and treatment for transcriptome analysis.

RNA sequencing
For transcriptome sequencing (RNA-seq), three biological replicates were collected from both the control and IBA-treated groups at each time point.In total, 24 RNAseq libraries (two treatments × four time points × three biological replicates) were generated.Total RNA was isolated using TRIzol reagent, and the libraries were constructed and sequenced on an Illumina HiSeq™ 2500 platform at BMK Company, Beijing, China.Raw sequence reads, comprising 150 bp paired-end reads, were filtered and aligned as previously described by Ahmad and colleagues [33].

Sequence alignment to the Morus notabilis genome and RNA-sequencing data analysis
Following sequencing, high-quality reads were obtained by removing adapter sequences, low-quality reads, and ambiguous nucleotides (N).Concurrently, during the trimming and filtering process, descriptive statistics for the resultant high-quality data were calculated, including Q20, Q30 scores, GC content, and the level of sequence duplication.These high-quality reads were then used for further analysis.They were aligned to the Morus notabilis reference genome available at the Morus notabilis reference genome (https://morus.biodb.org/browse)using HISAT2 software [34].DEGs were identified by calculating the log2 fold-change (FC) of gene expression at different treatment stages.DEGs were selected based on |log2(FC)| ≥ 2 and a statistical significance threshold of P ≤ 0.05.The DESeq tool in R was employed to detect DEGs using the criteria of |log2 ratio| ≥ 1 and an adjusted P-value (false discovery rate, FDR) ≤ 0.05 [35,36].
The KEGG (http://www.kegg.jp,accessed on 13 October 2023) and GO (http://geneontology.org, accessed on 28 October 2023) databases were utilized to perform enrichment analyses of transcripts and DEGs per sample.KEGG facilitates the prediction of protein interaction networks and their functions in various cellular processes.GO enrichment analysis was applied to categorize the primary biological functions of the DEGs in terms of molecular functions, cellular components, and biological processes.The hypergeometric test was used to identify significantly enriched pathways and GO terms among the DEGs compared to the genomic background.The resultant P-values were adjusted to control the FDR, with an FDR ≤ 0.05 considered significant.The DESeq R package was employed to apply the hypergeometric test for enrichment analysis.

Weighted gene co-expression network analysis
WGCNA was conducted based on the expression correlation patterns among DEGs.The DEGs were analyzed using the log2-transformed FPKM values plus one as input, and the soft thresholding power was determined by the scale-free network criterion [37].The lowest power value at which the scale independence reached a plateau (or exceeded 0.8) was chosen for downstream analysis, and the changes in gene connectivity at various power values were also examined [38][39][40].Genes were clustered into modules using dynamic tree cutting.A gene clustering dendrogram was constructed based on gene expression correlations, and gene modules were defined according to the clustering dendrogram.Modules with similar expression profiles were then merged based on the similarity of their module eigengenes, with a minimum of 50 genes per module and a merging threshold of 0.8.Modules were identified as significant through module eigengene analysis, and relevant modules were selected for more detailed investigation.

Validation of DEGs by RT-qPCR
Real-time quantitative PCR (RT-qPCR) was used to validate the transcriptome data, and 10 DEGs were randomly selected for RT-qPCR.The Actin gene was used as an internal reference gene [41].Specific primers: 5 Primers were designed according to the sequences of each gene (Table 1), and the differential gene expression levels were detected on the Bio-Rad Real-Time Fluorescence Quantitative PCR Instrument.The RT-qPCR system consisted of 2 × RealStar Green Fast Mixture 10 µL(Genstar, Beijing, China), template cDNA 1 µL, forward/reverse primers 0.5 µL (10 µmol/L) each, ddH2O 8 µL PCR program: 94 ℃ predenaturation 2 min; 94 ℃ denaturation 15 s, 60 ℃ annealing30 s, cycling 40 times.The relative gene expression was calculated by the 2 ΔΔCt method [42].

Sequencing data quality control and comparison of reference genomes
The 24 RNA-seq libraries that were generated underwent analysis, and the short sequences produced through sequencing constituted the raw data.Given that RNA extraction, library preparation, and sequencing can introduce redundant or low-quality data, clean data were acquired by filtering the raw data to remove duplicated reads, reads containing adapters, reads with a high proportion of N, and low-quality reads.This process involved quality assessment and control.From the 24 samples, a total of 150.20 Gb of valid data was obtained, averaging 5.02 Gb per sample, with an individual sample data size of around 7.14 Gb on average.These data were saved in the FastQ file format, facilitating the smooth progression of subsequent bioinformatics analyses.The data output statistics for each sample are presented in the accompanying table (Table 2).Post sequencing quality control, a  3).

Repeat relevance assessment
To identify differentially expressed genes of genuine interest, it is necessary to account for and mitigate the impact of this biological variability.In this study, the correlation between gene expression levels across samples serves as a critical metric for evaluating the reproducibility of the biological experiments, confirming the validity of the identified differentially expressed genes, and aiding in the identification of outlier samples.We employed Pearson's correlation coefficient (r) as the measure of correlation between biological replicates [43], with an r2 value approaching 1 denoting a strong correlation between two replicate samples.The heatmap depicting the sample correlations in this study is presented in Fig. 2. The results indicated that most replicates clustered together, suggesting that the biological replicates were generally well-established.The highest correlation was observed between the treatment groups across different time points, while a lower correlation was apparent within the CK, which also suggests that the gene expression in the samples underwent significant changes posttreatment, consistent with the expected pattern of the experiment.

Analysis of EDGs in treatment and expression groups
Before and after treatment with 1000 mg/L groups, suggesting that these DEGs continued to exhibit significant changes throughout the course of the treatment with IBA.
In parallel, this study utilized STEM software to normalize the expression data based on log2(fpkm + 1) using the expression counts.This normalization was then used to analyze the expression trends of the 4,124 DEGs in both the CK group and the treatment group, Fig. 4 shows two different trends The analysis revealed two distinct expression patterns with significant differences in expression levels among these DEGs in the treatment group, while five distinct expression patterns with significant differences were observed in the CK group.The trends in expression of these DEGs in the IBA treatment group and the CK group were markedly different.

GO and KEGG analysis of DEGs
To elucidate the molecular mechanisms in mulberry following IBA treatment, GO and KEGG analyses were conducted on differentially expressed genes (DEGs) from four temporal stages.GO enrichment analysis of the top 20 categories (Fig. 5) categorizes DEGs into Diverging from the GO annotations, the KEGG pathway enrichment results varied across the four stages (Fig. 6).In the comparison between CK-A and IBA-1, the pathways of photosynthesis -antenna proteins (ko00196), taurine and hypotaurine metabolism (ko00430), and caffeine metabolism (ko00232) were significantly enriched.In the comparison between CK-B and IBA-2, pathways including biosynthesis of various secondary metabolites (ko00998), flavonoid biosynthesis (ko00941), and taurine and hypotaurine metabolism (ko00430) were notably enriched.For CK-C vs. IBA-3, pathways such as glycosphingolipid biosynthesis -lacto and neolacto series (ko00601), taurine and hypotaurine metabolism (ko00430), and flavonoid biosynthesis (ko00941) were significantly enriched.Finally, in the comparison between CK-D and IBA-4, pathways like taurine and hypotaurine metabolism (ko00430), zeatin biosynthesis (ko00908), and synthesis and degradation of ketone bodies (ko00072) were significantly enriched.Collectively, these results suggest that most DEGs were associated with secondary metabolite synthesis pathways under IBA treatment, indicating a strong engagement in organic matter metabolism.Remarkably, after analyzing the number of genes enriched in each pathway, the pathway with the highest number of enriched genes was plant hormone signal transduction (Table 4), further suggesting that the phytohormone pathway is significantly activated following IBA treatment.This observation aligns with the findings from the aforementioned ARs production study and

Analysis of differential expression related to hormone signaling pathways in DEGs
As indicated previously, the phytohormone signaling pathways were markedly enriched in mulberry following treatment with 1000 mg/L − 1 IBA.Subsequent analysis revealed that differentially expressed genes (DEGs) were predominantly enriched in auxin (Fig. 7A), gibberellin (Fig. 7B), ethylene (Fig. 7C), brassinosteroid (Fig. 7D), and salicylic acid (Fig. 7E) pathways.The findings demonstrate that these hormone signaling pathways were active in the initial stage, with the majority of genes within these pathways being upregulated during this early phase.These genes exhibited significant upregulation in the initial period, leading to the conclusion that the transition from the commencement of treatment to the tissue healing stage is critical under the influence of 1000 mg/L − 1 IBA.Notably, in the auxin signaling pathway, AUX/IAA gene expression was highly active, with numerous genes showing an upregulation of more than eightfold in the initial period.In contrast, the TIR1 gene did not exhibit significant upregulation until the fourth period, which may be attributed to the accumulation of TIR1 in the plant CK group.In the gibberellin signaling pathway, the pattern was distinct from the other hormone pathways; the gibberellin receptor GID1, DELLA proteins, and TF were all downregulated during the initial period, with expression levels decreasing by 1-2-fold.This suggests that DELLA proteins act as negative regulators in the gibberellin signaling pathway, thereby inhibiting plant growth and development.

Transcription factor analysis of EDGs
The analysis of transcription factors among DEGs from the four periods, using Arabidopsis as the reference species, yielded a total of 410 annotated DEGs, with the top 20 by number showcased in Fig. 8.The RLK-Pelle_ DLSV category contained the highest number with 54 members, followed by AP2/ERF-ERF with 39 members, and MYB with 33 members.The categories hsf, RLK-Pelle_CrRLK1L-1, and RLK-Pelle_L-LEC had the fewest members, each with only 11.It was confirmed that most of these transcription factors are closely related to plant hormones.These findings corroborate the earlier results and collectively reinforce the connection between these DEGs and plant hormones.

Identification of key transcription factors of DEGs by WGCNA Soft threshold determination and clustering of genes in gene coexpression networks
From Fig. 9A, it is evident that a soft threshold value β of 22 resulted in a scale-free network fitting index R2 greater than 0.8, and the mean connectivity approached zero.This suggests that employing a β value of 22 enables the generation of a scale-free network that satisfies the analytical criteria; hence, β = 22 was selected for the construction of the scale-free network.
A dendrogram was generated based on the pairwise correlation of gene expression profiles (Fig. 9B).The dendrogram was truncated using dynamic tree cutting, grouping genes with similar expression patterns into the same branches, with each branch representing a distinct coexpression module.Following the amalgamation of modules with analogous expression patterns based on a threshold module similarity of 0.8, eight coexpression modules were delineated.Each module is denoted by a unique color, and the genes not classifiable into any module are represented by the color gray.The black module comprises the largest number of genes, totaling 1,081, succeeded by the blue module with 438 genes, the brown module with 192 genes, the green module with 179 genes, and the pink module with 133 genes.The magenta module includes 116 genes, the purple module comprises 105 genes, and the green-yellow module contains the smallest number of genes, with just 80 genes.The obtained modules were analyzed for correlation with the samples, and eight modules related to different varieties and treatment times were obtained.Some of the modules were highly correlated with treatments and periods (Fig. 9C).By observing the correlation between the modules and the samples, it was found that magenta, black and brown modules were significantly positively correlated with the traits.Specifically, black and brown modules exhibited the strongest correlation, so the two modules, black and brown, were taken as the IBA treatment-related specific modules for in-depth analysis to excavate the core genes in the modules.

Promoter analysis of ARs development related transcription factor genes
To further investigate the relationship between transcription factors and genes involved in plant hormone biosynthesis and signal transduction, we performed cis-acting regulatory element prediction analysis on the sequences approximately 2000 bp upstream of these genes using the PlantCARE database, the results are shown in Fig. 12.As expected, in the black and brown modules, nearly all co-expressed hormone-related genes and promoters of genes associated with ARs formation contained cis-elements such as ERE-motif, MYB-motif, G-box, W-box, MYC-motif, and ARR-motif.The repeatability of these elements suggests that these genes are not just acting individually, but are more likely to be closely related to the regulation of root growth.

Real-time PCR validation of DEGs
To validate the RNA-seq results, 10 DEGs were randomly selected for RT-qPCR verification.The expression patterns of these genes were consistent between RNA-seq and RT-qPCR analyses, the result is shown in Fig. 13, which confirmed the accuracy and scientific validity of our experiment.

The expression of endogenous hormones in plants is closely related to the formation and development of ARs
Cutting is a widely used method of vegetative propagation in horticulture.The success of this method hinges on the regenerative capacity of plant tissues.Plant hormones play an essential regulatory role by interacting with transcription factors and other regulatory elements to direct cell division, morphogenesis, and functional differentiation, ultimately leading to the regeneration of roots and shoots following callus formation [44].In this study, a structure resembling a callus formed at the wound site of the cuttings, from which ARs subsequently developed.Previous research has indicated that endogenous hormones do not act independently in cuttings; rather, a synergistic action of multiple endogenous hormones is necessary to promote rooting and the initiation of shoots [45].Among these, IAA is a pivotal regulator of ARs formation post-cutting, while cytokinins, jasmonic acid, gibberellins, brassinosteroids, ethylene, and other hormones are known to induce or enhance the initiation of root primordia and the formation of ARs through interactions with auxins and cytokinins [46].
In the hormone-related gene expression analysis of mulberry cuttings, differentially expressed genes were identified across the five hormonal pathways-auxin, gibberellin, ethylene, brassinosteroids, and salicylic acid.This finding suggests that these hormones participate to varying extents and through diverse mechanisms in this complex biological process.The endogenous hormoneassociated genes related to IAA, gibberellin, and BR were most prevalent, with their biosynthesis and signal transduction genes upregulated across all four cutting time points.This upregulation suggests an increase in their levels and signaling activity, which may trigger cellular differentiation in mulberry leaves and the formation of ARs in cuttings.Previous research has shown that the balance between growth-promoting and growth-inhibiting hormones influences the formation of ARs [47].In this study, the biosynthesis genes for IAA and gibberellin were found to be upregulated at day 10 post-cutting, with a greater number and a higher ratio of upregulated genes for IAA than for gibberellin, implying that a high IAA/ gibberellin ratio may be favorable for ARs development in mulberry.Additionally, BR signaling genes primarily facilitated the formation of ARs in mulberry cuttings through signal transduction.The expression patterns of the BR signaling gene BRI1 and the gibberellin biosynthesis gene KAO2 were consistent in the later stages of cutting, suggesting that BR may enhance gibberellin biosynthesis by upregulating BRI1 expression.This aligns with findings that BR induces the expression of GA20ox2 in rice seedlings by increasing the levels of the signal transduction protein BRI1, which raises the levels of bioactive gibberellin and promotes the formation of ARs [48].

Transcription factors WRKY, NAC, MYB and TCP are closely related to the formation and development of ARs
Studies have shown that WRKY75, a member of the WRKY gene family, regulates the activity of phosphatases at the transcriptional level, thereby influencing the dynamics of auxin transport and lateral root development [49].In Arabidopsis thaliana, WRKY23 has been identified as promoting localized flavonol biosynthesis, which facilitates root growth and maturation.This process is regulated by auxin through the transcriptional responses of ARF7 and ARF19 [50].The wheat transcription factor TaNAC2-5 A binds directly to the promoter regions of genes encoding nitrate transporter and glutamine synthetase, enhancing root proliferation and the rate of nitrate uptake.This, in turn, improves nitrogen acquisition and ultimately increases grain yield [51].In Arabidopsis, the transcription factor NAC1 is specifically induced by wounding in leaf explants and aids in the regeneration of the root apex [52].In the same species, the MYB family gene MYB77 potentially affects the number of lateral roots through its interaction with ARF7 [53].Furthermore, MYB15 expression is initially upregulated in primary callus and later downregulated during rerooting, indicating its role in regulating callusinduced differentiation in tea plants [54].Aguilar-Martinez et al. [55] reported that AtBRC1/TCP18 interacts with both auxin and strigolactone pathways, contributing to the regulation of plant branching.In cotton, the type I protein GbTCP, analogous in function to AtTCP15, has been shown to reduce jasmonic acid levels and inhibit fiber elongation when silenced, while its overexpression in Arabidopsis thaliana promotes root hair elongation [56].Notably, in mulberry, these transcription factors are significantly associated with auxin within the network, suggesting a primary interaction with relevant genes during the induction and elongation phases of ARs formation in tea cuttings and playing a role in the regulation of ARs regeneration.Promoter analysis supports this hypothesis, indicating a necessity for further experiments to elucidate the molecular mechanisms and physiological roles of these regulatory factors in ARs regeneration.

Plant hormones work together to induce changes in plant roots
An exploratory analysis of the cis-regulatory elements associated with 18 transcription factors revealed that in the black and brown modules, nearly all co-expressed hormone-related genes and promoters of genes  [57,58].Notably, ABRE-motifs, involved in the abscisic acid response, were identified in gene18390 and gene19492, indicating a strong association with abscisic acid [59].The AuxRE motif, essential for the auxin pathway, was found in gene17435, with auxin response factors selectively binding to the AuxRE motif to regulate auxin signals [60].These findings suggest that in addition to governing downstream genes of this pathway, the endogenous hormones associated with ARs may participate in potential cross-regulation, possibly influencing the expression of genes within the plant hormone signal transduction pathway.

Conclusion
By analyzing RNA sequencing data from mulberry phloem before and after IBA treatment, we identified differentially expressed genes.Subsequently, we discovered 18 transcription factors that potentially promote ARs formation.These transcription factors, belonging to the WRKY, NAC, MYB, and TCP transcription factor families, are closely associated with the development of ARs.Exploratory analysis of the cis-regulatory elements associated with these transcription factors suggests their potential joint regulation of plant hormone signal transduction and a ARs formation.These families may play a pivotal role in the generation of ARs, warranting further study and exploration.

Fig. 1
Fig.1Changes in root morphology of mulberry across four periods

Fig. 4 Fig. 3 Fig. 6
Fig. 4 Trend analysis of gene coexpression of all DEGs over eight periods in two treatments.All DEGs were categorized into two expression trends (A, B)

Fig. 5
Fig. 5 GO enrichment analysis of DEGs from four time periods

Fig. 7
Fig. 7 Expression of DEGs in plant hormone pathways

Fig. 8
Fig. 8 Transcription factor analysis of all DEGs in four time periods

Fig. 9
Fig. 9 Presents a WGCNA of the gene expression matrix in mulberry.(A) The most appropriate soft threshold was determined by plotting scale independence and mean connectivity.(B) A dendrogram based on coexpression network analysis depicts the hierarchical clustering of genes, with the module colors represented on the X-axis.(C) The module-sample association is shown, where each row corresponds to a module color-coordinated with that in part B, and each column represents a sample.The correlations between the various modules are indicated by the values inside the colored boxes

Fig. 10
Fig.10 The top 150 connectivity gene networks in the black module.The orange color indicates the key transcriptional genes that have been screened

Fig. 11 Fig. 13 Fig. 12
Fig.11The top 150 connectivity gene networks in the brown module.The orange color indicates the key transcriptional genes that have been screened

Table 2
Statistical tables of sequencing data

Table 3
Statistics of sequence comparison results of sample sequencing data with selected reference genomes regulated DEGs were obtained; CK-D vs. IBA-4 identified 901 DEGs, of which 572 were upregulated and 329 were downregulated.The comparisons CK-A vs. IBA-1 and CK-B vs. IBA-2 shared 393 DEGs, CK-B vs. IBA-2 and CK-C vs. IBA-3 shared 124 DEGs, CK-C vs. IBA-3 and CK-D vs. IBA-4 shared 78 DEGs, and CK-A vs. IBA-1 and CK-D vs. IBA-4 had 341 DEGs in common (Fig. 3).There were six genes common to all four comparison

Table 4
Statistics of the number of DEGs enriched in the KEGG pathway

Table 5
Functional annotation of core transcription factors in correlation-specific modules